Visibility of Cold Atomic Gases in Optical Lattices for Finite Temperatures 
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In nearly all experiments with ultracold atoms time-of-flight pictures are the only data available. 
In this paper we present an analytical strong-coupling calculation for those time-of-flight pictures 
of bosons in an optical lattice in the Mott phase. This allows us to determine the visibility, which 
QQ , quantifies the contrast of peaks in the time-of-flight pictures, and we suggest how to use it as a 

. thermometer. 
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I. INTRODUCTION 



Systems of ultracold bosonic gases in optical lattices have recently become a popular research topic [1, 2]. After 
their theoretical suggestion [3, 4] and first experimental realization [5] it soon became clear that they represent model 
systems for solid-state physics with a yet unprecedented level of control. Both the periodic one-particle potential, which 
. is superimposed with an additional harmonic trap for the purpose of confinement, and the two-particle interaction 
strength can be experimentally tuned with high precision [6] . With this it is now possible to achieve strong correlations 
in these systems even though their particle densities are more than five orders of magnitude less than for air. In 
addition, due to the absence of impurities, they are viewed as idealized condensed matter systems which allow for a 
clear theoretical analysis [7, 8]. 
■ Current research in optical lattices is driven by a number of cornerstone experiments which investigate, among 



other things, lower dimensionality [9], the formation of a distinct shell structure in the particle density [10], and a 
more quantitative analysis of the time-of-flight interference patterns via their visibility [11]. Furthermore, dynamical 
aspects are now under crucial investigation [12] as, for instance, the collapse and revival of the interference patters 
, in nonequilibrium situations [13-15]. Furthermore, using either bosonic or fermionic atoms in optical lattices reveals 
(~| a contrasting bunching or antibunching behavior which can be fully attributed to the different quantum statistics of 
O '. each atomic species [16]. Mixtures of bosons and fermions are also currently studied [17-19], and this has even led 
^ to the creation of heteronuclear molecules [20]. Thus, in the near future it will be possible to test the theoretical 

[ prediction of new phases in optical lattices arising from an additional dipole-dipole interaction [21]. Another emerging 
^— I I line of research is the investigation of the effect of controlled disorder which can be created by several methods [1]. It 
^ ' turns out that an additional frozen random potential leads to a Bose-glass phase due to the localization of bosons in 
, the randomly distributed minima [3, 22-24]. Finally, it has been proposed to use ultracold bosons in optical lattices 

■ to realize a quantum computer [25-27]. 

^— ^ , An optical lattice is created by pairs of counter-propagating laser beams in all three dimensions. Loading bosons 
with mass M in this laser field, the Stark effect leads to an effective periodic one-particle potential V^xtCx) with 
[ lattice spacing a — X/2 and laser wave length A. Thus, the energy scale of the system is set by the recoil energy 

■ Efi = ir'^h^ /2Mo? , which is used to introduce dimensionless energies E = E/Efj. Neglecting the overall harmonic 
I trapping potential, which is superimposed in order to spatially confine the system, spin-polarized bosons in an optical 

• • ■ lattice can be described within the grand-canonical ensemble by the Bose- Hubbard Hamiltonian [3, 4] 



where aj and a.; denote the standard creation and annihilation operators at site i, /i denotes the chemical potential, 
and the sum over < i, j > includes only pairs of nearest neighbors. Both the hopping matrix element J and the on-site 
interaction U define the major energy scales of the system, which turn out to depend crucially on the strength Vq of 
the laser field (see Appendix A). If Vq is varied, the bosons can undergo a quantum phase transition. For small Vq the 
hopping matrix element J is large and the bosons can tunnel from site to site to explore the whole lattice. This leads 
to a superfluid state which is characterized by long-range correlations, a continuous excitation spectrum, and a finite 
compressibility. In the situation of large Vq the hopping matrix element becomes negligibly small and the bosons can 
no longer tunnel to the neighboring sites, so the occupation number of the sites is fixed. This so-called Mott phase 
has no long-range correlation, shows a gap in the excitation spectrum, and is nearly incompressible. In time-of-flight 
pictures these two phases can also be distinguished by their distinct interference patterns. In the superfluid phase the 
atoms are delocalized over the complete lattice. Therefore, according to the Heisenberg principle, their momentum 
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uncertainty is small, which leads to sharp Bragg peaks in the time-of-fiight pictures. In the Mott phase, on the other 
hand, all particles are strongly localized at lattice sites, leading to a large momentum uncertainty and ultimately to 
a uniform cloud during the expansion. 

It is of particular interest to determine how the location of the transition from the superfluid to the Mott phase 
depends on the respective system parameters. Usually, one assumes that the temperature in the experiments is so 
low that thermal effects are completely negligible. In that case the phase boundary has been determined analytically 
within both a mean-field theory [3, 28] and a strong-coupling approach [29, 30] as well as numerically by Monte-Carlo 
simulations [31-33]. Recently, with the field-theoretic concept of effective potential, the complete quantum phase 
diagram has been calculated essentially exact in excellent agreement with available numerical data [34]. Only some 
theoretical work has been initiated to include thermal effects in a systematic way [22, 35-37]. However, until today, a 
reliable method of determining the temperature of the bosons in an optical lattice is not known [16]. Therefore, more 
experimental and theoretical studies are needed that aim at designing a thermometer for these systems. 
This motivates the present paper where we investigate how the temperature affects the time-of-flight pictures for such 
a lattice system. We start in Section II with determining perturbatively the correlation function for the homogeneous 
Bose-Hubbard model (1) within the Mott phase. To this end we work out a hopping expansion up to second order 
at finite temperatures. Note that this hopping expansion is closely related to the random walk expansion worked out 
in Refs. [38-40]. The resulting correlation function is used in Section III to qualitatively reconstruct time-of-flight 
absorption pictures which are taken after switching off the one-particle potential. From this we calculate in Section IV 
the visibility, which quantifies the contrast of the time-of-flight pictures, determine how it changes with varying lattice 
depth Vq, and compare our results with experimental data. The Appendix summarizes more technical material. In 
Appendix A we determine how the Bose Hubbard parameters J and U depend on the lattice depth Vq. In particular, 
we show that the commonly used harmonic approximation [41, 42] deviates significantly from numerically determined 
results for large lattice depth. Finally, we work out in Appendix B how an additional external harmonic confining 
potential fixes the average number of bosons per site. 



II. CORRELATION FUNCTION 



The time-of-flight absorption pictures rely on the correlation function 

{alaj)=Z-'Tr{a]dje-0"} (2) 

with the partition function Z = Tr As it is not possible to calculate this quantity analytically by ex- 

actly diagonalizing the Bose-Hubbard Hamiltonian, we have to employ a perturbative scheme. Therefore, we 
will restrict ourselves to the strong- coupling regime of the Mott phase, where the exactly solvable on-site part 

Hq = [■j?7,,(rii — 1) — fifiij, with the number operator r)^ = a|dj, determines the unperturbed system and the 
hopping part V = —J^^ijy afoj can be considered a small perturbation. With this we can rewrite (2) using the 
imaginary-time evolution operator in the Dirac picture IJt>{t, t') = pHar /h^H{T -T)/h^-HoT jh according to 

(ata,) =Z-iTv{a|a,e-'5^°{7D(/i/3,0)} , (3) 
where also the partition function can be written in term of C/d: 

2 = Tr{e-'^^''Z7D(?i^,0)} . (4) 
The Dirac time-evolution operator can now be expressed by the Dyson series: 



?7D(r,r') = 1 + ^ y dTiyD(Ti) + J dn J dT2Vn{n)VB{r2) + ■ 



(5) 



With this we obtain from Eq. (4) an expansion of the partition function in powers of the tunnel matrix element J. 
As the trace is only sensitive to the diagonal part and as one always needs in a cubic lattice an even number of steps 
to return to the starting point, all contributions with an odd power of J must vanish leaving 

Z = ^(0) + J2Z(2) + . . . (6) 
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with the coefficients 



zW=Tr{e-'5^''} , 

fhf3 Ti 
e-""'' J dn J dT2Mn)Mr2) 



(7) 
(8) 







The correlation function (3) is now calculated in the same manner by applying the Dyson series (5), so that we end 
up with a perturbation series in J: 



Here the respective expansion coefficients read up to the second order in J: 
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Tr I alcne-^"' 



ft/3 



(n) 



.(2) _ 
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Tr < alaje 



/f c 
dn J dr2Vb(Tl)Vb(T2) \ - 



(9) 

(10) 
(11) 

(12) 







Now we have to explicitly evaluate all the respective traces. To this end we use the occupation number basis |n) = 
Hi l^^i), in which the unperturbed system is diagonal, i.e. Ho\n) — Eftlfi) with Eft = with the unperturbed 

single-site energies of the system given by En = Un{n — l)/2 — \xn. Thus, representing the trace according to 



n 

the coefficients of the partition function (7) and (8) reduce to 

^(0) ^^g-/3E,i5n, ^ 



z(^) =/3 E E 



ni{n^ + 1) 



<ij> n 



En,+1 + En--l — En- — 



(13) 

(14) 
(15) 



Analogous calculations yield the coefficients of the correlation function (10)-(12). With the notation 



y 



these can be written as 

cl? =k,3{ni)i 



(1) / 2ni{ni + 1) 

U \{ni-nj + l){nj-ni + l)/^ 



.(2) 



^d{i,j),2 + 2(5d(j,j),\/2 / Uiirij + l){ni + 1) 



n^rijini + 1) 



+ 



U 



{uj -ni + l)(n; - m + l) {m - n, + l){ni - Uj + ^) / 

nf{ni + 1) \ / \ / '^ii'^i + 1) 



nj - rij + 1 / . , 



(16) 

(17) 
(18) 

(19) 



where d{i,j) is the distance in the Euclidean norm measured in units of the lattice spacing and = 6 is the number 
of next neighbors in a three-dimensional cubic lattice. 
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III. TIME-OF-FLIGHT ABSORPTION PICTURES 



Measurements in most experiments on ultracold atomic gases are made as follows. The trapping potential is switched 
off allowing the gas to freely expand during a short time of flight t. Then an absorption picture is taken, which maps 
the particle density in real space to a plane. Due to the diluteness we can neglect any interaction between the atoms, 
so the particles will approximately move with a constant velocity given by their momentum at the moment of release 
from the trap. Furthermore, we will assume that the condensate in the trap is point-like, so that the distribution in 
real space is a mirror of the distribution in momentum space connected by the relation Kk. = mr/t. The momentum 
space distribution in the optical lattice is given by [43] 



n(k) = mk)|25(k) , 

where the quasi-momentum distribution reads 



with the Fourier transform of the Wannier function being defined according to 

dh 



With (3) and (17)-(19) we get the quasi-momentum distribution S'(k) as 



5(k, T) = So{T) + 2^Si{T) J2 cos(/cia) + 2^S2{T) 



i=l 



-3 -f ^ 2cos(A;ia) cos(A;ja) 



+ 



with the temperature dependent coefficients 



5o(r) =Ns{nk)k + Ns- 



U 



nk 



nkjni + 1) 
m-nk + l/kj 



{nk)k 



nkjni + 1) 
ni-Hk + l 



k,l 



Si{T) =Ns 
S2{T) =Ns 



2{nk + l)nfe 



(rife -ni + l)(n/ - Ufe -I- 1) / 
nfc(nz -h l)(nm + 1) ^ 



nkniirim + 1) 



{m ~ nk + l){nm -Uk + 1) {tim -Uk + l){nm ~ni + 1) / 



(20) 
(21) 

(22) 



(23) 



(24) 
(25) 
(26) 



Here Ns denotes the number of lattice sites. Taking an absorption photograph projects the density distribution onto 
a plane. Theoretically, this corresponds to integrating the particle density over the 2-axis 



n{x,y,t) = 



— I I dzn 
ht J 



/Mr 
\lit 



(27) 



Together with a factorization of the Wannier function w(k) = nj=i '^{^i) this yields the following result for the 
time-of-flight pictures: 



n{x,y,t) 



2Max 



w 



cos ■ 



-+COS 



\ ht J 

2May 



V ht 



So{T) + 2^S^iT) 



Max May , 
cos h cos h W{a) 



ht 



ht 



2^52(T) 



ht ht 

where we have introduced the abbreviation 



r,.,r. X . Max May , Max^^^, , , May„^, , 
W{2a)+A cos — — cos + 4 cos ——W(a) + 4 cos —-^W(a] 
^ ^ ht ht ht ^ ' ht ^ ' 



W{1) 



oo 

_ M r 

~ht J 



w 



fMz\ 



cos 



Mh 



dz . 



■+... ,(28) 



(29) 
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FIG. 1: Time-of-flight pictures for Vo = 8 (column a), Vb ~ 14 (column b), Vo = 18 (column c), and Vo — 30 (column d). 
[second] row shows first- [second-] order calculation and the last row corresponds to the experimental data of Ref. [44]. 



First 



Using the Wannicr function (A4) of the optical lattice within the harmonic approximation yields the Fourier transform 
[11]: 



\Mh)\' 



= exp 




This allows to calculate (29) explicitly 



(30) 



W{1) = e 

In the following we present theoretical density plots of (28) only for the T — > limit, i.e. 



^o(O) 
^i(O) 
52(0) 



Nsn, 

2Nsn{n + 1) , 
Nsn{n + l)(2n - 



1): 



(31) 



(32) 
(33) 
(34) 



as a procedure to determine the true temperature of bosons in an optical lattice experiment is not yet known. Fig. 1 
compares this result to the experimental findings in Refs. [11, 44]. There the laser had the wavelength A ~ 850 nm 
and the waist w = 130 /im in order to produce a cubic optical lattice, which was superimposed with an external 
harmonic potential of frequency Um — 2tt x 15 Hz. This configuration was filled with N — 2.2 x 10^ Rubidium-87 
atoms with an s-wave scattering length of qbb = 5.34 nm, so the recoil energy is = 2.10 x 10^"^" J. In Appendix 
B we determine that this experimental situation corresponds to the occupation number of n = 2 bosons per site. 
For the high laser intensities in columns (c) and (d) of Fig. 1, where the perturbation parameter J/U is sufficiently 
small, the theoretical pictures match the experimental ones, and there is nearly no difference between the first and 
the second hopping order. In columns (a) and (b), on the other hand, the theory does not fit the experiment. In the 
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theoretical picture we can observe the formation of a central and various neighboring peaks, which become sharper 
in the second order, but the neighboring peaks especially are very faint and, furthermore, some additional unphysical 
peaks appear. These problems arise because the perturbation parameter J /U is no longer a small quantity in the 
superfluid regime, so that our approximation is no longer valid. Furthermore, it becomes necessary to introduce the 
condensate as an order parameter and to take into account its influence on the quasi-momentum distribution. These 
deficiencies can be removed by, for instance, calculating the efii'ective action within the hopping expansion [45]. 



IV. VISIBILITY 



Following the approach of Refs. [11, 44], the contrast of peaks in the time-of-flight pictures is quantified by the 
visibility, which is defined as 

, ^max ''^min /or\ 
V = . (35) 

Here one considers the maximum rimax at the first side peak of the integrated density (28), while the minimum rimin 
is taken with the same distance to the central peak. 

We start with observing that, using the harmonic approximation, the Wannicr envelope (30) cancels out in (35), so 
that the result only depends on the quasi-momentum distribution (23). We find both the maximum and the minimum 
by evaluating (23) at the points x = 2Trht/Ma,y = and x = y = ^Trht/Ma, respectively. With this both the 
numerator and denominator of the visibility (35) yield in second order of the ratio J/U : 



^ ^ 2(1 - cos ^/27^)5l(T)^ + 8 sin^ V2TrS2{T)^ + ... 

5o(T) + 2(1 + cos \/27r)5i(T)^ + 4(1 + cos \/27r)S2(T)^ + . . . ■ 

Expanding this fraction in powers of zJ/U yields 
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with the prefactors 



V = v^{Tf-^+V2iT)r-l] +... (37) 



.r(T) = '-^^^m , (38) 



sm^ v^TT 

V2{T) = 



S^iT) (S,{T) 



So{T) \So{T) 



(39) 



which reduce at T = to 



.,(0) = l^^^|(^(n + l)>0, (40) 

^2(0)= -^";'^" (n+l)<0. (41) 

With this we have generalized Eq. (5) from Ref. [44] to finite temperatures with a corrected prefactor. Before we 
discuss the impact of temperature on the visibility, we compare in the T — > limit Eq. (37) with direct calculations 
from Eq. (28) in Fig. 2. In the double logarithmic plot we use again the occupation number of n = 2 bosons per site. 
All curves show a linear behavior in the Mott phase, while in the superfluid regime our direct calculations produce 
the unphysical result of the visibility being larger than 1. This is, again, due to fact that the above approximations 
are not appropriate in this parameter region. In the Mott region, on the other hand, all curves coincide and the linear 
dependence in zJ/U is supported by experimental data [11, 44]. However, we note that the theoretical prefactor 
(40) turns out to be too small by a factor of about 2. This discrepancy may arise from the approximation that our 
calculation is carried out for the homogeneous case, so the influence of the shell structure on the visibility is neglected. 
However, taking this effect into account should decrease the theoretical result even further, as the correlation function 
depends in a nonlinear way on the occupation number. In addition, we have also neglected the existence of superfluid 
regions between the Mott shells, but these are small and highly depleted deep in the Mott regime, so this effect should 
not play a decisive role. 
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FIG. 2: Comparison of the visibility in first order (dashed) and second order (solid) using Eq. (36) to the first-order expanded 
result (dotted) given in Eq. (37) for T = using the harmonic approximation for the Wannier functions. 
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FIG. 3: Comparison of the visibility in the harmonic approximation (dashed) and with numerical determined Wannier function 
in first order (solid) for T = to the experimental data (dots). 



Finally, we used the harmonic approximation for the Fourier transform of the Wannier functions so they cancelled in 
(36). As we show in Appendix A that the harmonic approximation can lead to large deviations, we investigate now its 
influence on the visibility. We read off from Fig. 3, that the curve with the numerical Wannier envelope in the Mott 
regime no longer has a linear dependence on J /U and lies above the experimental data. This discrepancy between 
our theoretical prediction for the visibility and the experimental findings can be explained as a temperature effect. 
To this end we investigate at first the expanded result with the harmonic approximation (37), where it is sufficient to 
study the prefactor (38) for different chemical potentials, which corresponds to different occupation numbers. We see 
in Fig. 4, that the thermal fluctuations slowly destroy the correlation and thus decrease the visibility. This qualitative 
finding within the harmonic approximation remains valid once the Wannier functions are calculated numerically. 
Furthermore, we observe in Fig. 3 an increasing discrepancy between our theoretical prediction and the experimental 
findings for larger lattice depths. Note that a detailed analysis shows that this tendency is also present in the 
experimental data [11]. This discrepancy suggests an increasing temperature of the bosons in the lattice. This would 
support the assumption that the ramping of the optical lattice leads to an adiabatic heating of the Bose gas [32, 46]. 
In order to take into account this adiabatic heating one should determine the relation between temperature T and 
the lattice depth Vq from the condition that the entropy remains constant throughout the ramping process. This task 
is certainly a nontrivial one as this adiabatic calibration curve T — T{Vo) should be determined across the quantum 
phase transition. We expect that such a calibration curve could be deduced in a consistent way within the above 
mentioned effective action approach [45] . Once the inhomogeneity of the trap is included, it would be possible to get 
a visibility curve for constant entropy, which could then be used to determine the entropy present in the experiment 
and, finally, to determine the temperature for certain system parameters. One could then use remaining deviations 
from the theory to check if any non-adiabatic heating or non-equilibrium final states occur. 



V. CONCLUSION AND OUTLOOK 



In this paper we calculated and analyzed the time-of-flight pictures and the resulting visibility for finite temperatures 
in a perturbative scheme for strong interactions going beyond the commonly used harmonic approximation. However, 
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FIG. 4: Prefactor given by Eq. (38) for n/U = 1.5 (solid), jji/U = 2.5 (dashed), and ji/U = 3.5 (dotted). 



we did not achieve a full understanding of the experimental data, as the experimental situation is quite more complex 
than initially presumed. Comparing our theoretical calculations with experimental data suggested that the adiabatic 
heating process during the ramping of the optical lattice has to be taken into account. Thus, it becomes necessary 
to find a calibration curve, which allows to determine the temperature T as a function of the potential depth Vq for 
a given entropy. 
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APPENDIX A: PARAMETERS OF THE BOSE-HUBBARD MODEL 



We start from the second-quantized Hamiltonian [4] 



H= j d^x^^x) 



^(X) + ^ / ^3^1 / d3x2V'^(Xi)V^t(x2)yi„,(xi,X2)V'(Xl)V'(x2) , (Al) 



where 0(x). V'^(x) are the usual bosonic field operators. The external potential is given by the optical lattice 
Vc^t{x.) = Vq sin^ (nxj/a), while we assume that the interaction is of the contact type Vint(xi,X2) = g 6{xi — X2) 

with the strength g = AnaBBfi^ /M , where obb is the s-wave scattering length. 

As the external potential T4xt(x) is periodic, we can expand the field operators in Wannier states w(x— Xj), which form 

a complete and orthonormal basis of functions localized at the respective sites i. As we are at very low temperatures, 
we restrict ourselves only to the lowest energy band. This is implemented by using the decomposition -(/"(x) = 
ai«;(x — Xj), which yields (1) with the matrix elements: 



U{i) 
J{i,j) 



47raBB?i^ 



j d^x\w{x- 
■ J d^xw*{x — Xj) 



t(x) 



W{x-Xj) 



(A2) 
(A3) 



where the coefficients U = U{i) and J = J{i,j) with i and j being next neighbors turn out to be independent of the 
sites, and the chemical potential is given by /i = /i' + J{i,i). Note that due to the orthonormality of the Wannier 
functions J is independent of the dimension while U is not, as the Wannier function factorizes in the respective spatial 
dimensions: 'w{x) = ni=i w{xi). 

If we use now the harmonic approximation [41, 42] and assume that the ground-state wave- function is given by the 
Gaussian 



w{xi) 




exp 



TT 



(A4) 
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FIG. 5: Comparison of the harmonic approximation (dashed) and the numerical results (plain) for the Wannier function (a), 
the on-site energy U (b), the hopping energy J (c) and the ratio of the last two (d) using the experimental parameter of the 
Bloch group. 



these integrals are directly solved, leading to 

J=(^!^^l)voe-^^/^/^ (A5) 

Alternatively, we can determine the Wannier functions numerically [47] and use this result to calculate the on-site 
interaction strength U via (A2). In contrast, the tunneling parameter J should not be determined from (A3) as 
numerical differentiation leads to low precision. Instead we use the fact that J does not depend on the dimension and 
the definition of the Wannier functions ^^{x — Xi) = Ng^^"^ Yl,k ^~^''^'^k{x), with the single particle Schrodinger wave 
function of the lowest band (f)k{x) and the corresponding eigenenergy E{k). With this (A3) reduces to [48] 

J = ^ e''''^''*-''i^E{k) . (A7) 

When we compare the analytical approaches with the numerical one in Fig. 5(a), we find, that the harmonic ap- 
proximation to the Wannier function nearly fits the numerical counterpart but for one special feature: a Gaussian is 
always positive while the numerical curve shows oscillations around zero. Note that it is, in fact, indispensable that 
the Wannier functions have also negative values in order to guarantee their orthogonality. This discrepancy leads to a 
significant deviation of the qualitative behavior in both Bose-Hubbard parameters over the whole energy scale. In the 
on-site energy U in Fig. 5(b) the harmonic approximation yields values, which are too large, but the relative error is 
at least slowly decreasing for higher laser intensities, while the relative error of the hopping parameter J in Fig. 5(c) 
is even growing. The ratio of these two given in Fig. 5(d) is the central calibration curve to connect the experimental 
data to the theory. Because of the higher precision, we use this numerical calibration curve throughout our paper. 
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FIG. 6: Shell structure slice of the Mott-phase for T = (solid), kaT/U = 0.05 (dashed), and kaT/U = 0.2 (dotted). 
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FIG. 7: Calibration curve for the chemical potential in the experiment of the Bloch group. 




APPENDIX B: SHELL STRUCTURE 



For our consideration it is important to know how many particles are at each site depending on the total particle 
number. In theory this number for T = is given by the next integer to the ratio l-i/U. 

In the homogeneous case, however, there is no natural restriction to the number of sites. Therefore, for any number 
per site one can find a number of sites, which gives the desired total particle number. In the experiments, however, 
we have an additional overall harmonic trapping potential with a frequency lo present. Additionally one has to take 
into account that the harmonic trap is enhanced by the shape of the laser beam. As the laser beam has a Gaussian 
intensity distribution, the lattice potential is deeper in the center of the laser beam than at the edge. The influence 
of this laser inhomogeneity can be described by an effective frequency [49] : 



8Vo - AEr^/VoJe^ 



(Bl) 



which limits the number of occupied sites. Additionally, it changes the local chemical potential and leads to concentric 

shells of Mott insulators, which have different numbers of particles per site. This is known as the wedding cake 
structure, which has been confirmed by comparing theoretical predictions [50] with experimental results [10]. A non- 
vanishing temperature blurs the boundaries between this shells as depicted in Fig. 6. As it is possible to measure 
this density variations one could use the shape of the boundaries to determine the temperature. But in this picture 
we neglected the superfluid shells, which are between the Mott shells. They are present even at T = and produce a 
similar blurring. 

If we neglect the discrete nature of the sites, determine the radii of the layers and subsequently add their volumes, 
we obtain for the equation of state [46] 



N : 




0<c<n/U 



(B2) 



This yields the calibration curve of Fig. 7, which has only negligible deviations from the exact sum over all lattice 
sites. From the calibration curve we can determine from a given particle number N the chemical potential ji for the 



11 



experiment and from this the occupation number in the central Mott shell. As the experiments of the Bloch group 
are done with a total particle numbers of about = 2.2 x 10^, the chemical potential is slightly smaller than 2U . 
Hence, the maximal occupation number is 2, which is consistent with the claims of Ref. [11]. 
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